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Abstract 

We have developed a cosmic ray (CR) shock code in one dimensional spherical 
geometry with which the particle distribution, the gas flow and their nonlinear 
interaction can be followed numerically in a frame comoving with an expanding 
shock. In order to accommodate a very wide dynamic range of diffusion length 
scales in the CR shock problem, we have incorporated subzone shock tracking and 
adaptive mesh refinement techniques. We find the spatial grid resolution required for 
numerical convergence is less stringent in this code compared to typical, fixed-grid 
Eulerian codes. The improved convergence behavior derives from maintaining the 
shock discontinuity inside the same grid zone in the comoving code. That feature 
improves numerical estimates of the compression rate experienced by CRs crossing 
the subshock compared to codes that allow the subshock to drift on the grid. Using 
this code with a Bohm-like diffusion model we have calculated the CR acceleration 
and the nonlinear feedback at supernova remnant shocks during the Sedov- Taylor 
stage. Similarly to plane-parallel shocks, with an adopted thermal leakage injection 
model, about 10~ 3 of the particles that pass through the shock and up to 60 % of 
the explosion energy are transferred to the CR component. These results are in good 
agreement with previous nonlinear spherical CR shock calculations of Berezhko and 
collaborators. 
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1 Introduction 



Collisionless shocks form ubiquitously in tenuous cosmic plasmas via collective 
interactions between particles and fields. The formation process of such shocks 
inevitably produces suprathermal particles, which can be further accelerated 
to become cosmic rays (CRs hereafter) through the interactions with reso- 
nantly scattering Alfven waves in the converging flows across a shock [11,28]. 
This so-called diffusive shock acceleration (DSA) is now widely accepted as 
the primary mechanism through which nonthermal particles are generated in 
a wide range of astrophysical environments. Both analytical calculations and 
numerical simulations have shown that DSA can be very efficient and that 
there are substantial and highly nonlinear back reactions from the CRs to 
the bulk flows and to the MHD wave turbulence mediating the CR diffusive 
transport (e.g. ,[28] and references therein). 

In the kinetic equation approach to numerical study of nonlinear CR acceler- 
ation at shocks, the diffusion-convection equation for the particle momentum 
distribution, f(p), is solved along with suitably modified gasdynamic equations 
[20]. A variety of analytic and semi- analytic methods have been developed 
to solve these equations for nonlinear, time-independent shocks (e.g., [26,9]). 
However, DSA modified shocks evolve as the CRs within them extend to higher 
energies over time. Many astrophysical shock systems evolve dynamically on 
time scales not particularly long compared to relevant DSA times, while the 
MHD waves that provide coupling between the CRs and the bulk plasma also 
should evolve (e.g. ,[32]). So, it is important to be able to model time depen- 
dent, nonlinear DSA in such systems. This will generally involve numerical 
simulations, a task that is quite challenging, partly because the full CR shock 
transition spans a very wide range of relevant length scales associated with 
the particle diffusion lengths, x^ip) = n(p)/u s , from CR injection scales near 
the shock to outer diffusion scales for the highest energy particles. Here, n(p) 
is the momentum dependent spatial diffusion coefficient and u s is the shock 
speed. The full shock structure forms on a scale determined by the highest 
momentum particles of dynamical significance; i.e. comparable to Xdipmax)- 
The greater of Xd(p max ) or the curvature radius of the shock defines the size 
of the system to be computed. On the other hand one must fully resolve flows 
on the scale of Xd{p m in) in order to follow properly the transport of the lowest 
momentum CRs. That scale will generally be only somewhat greater than the 
physical thickness of the viscous subshock that thermalizes most of the plasma 
flowing through the shock. These various scales may differ by many orders of 
magnitude for realistic models of n(p). 

Two approaches have so far been applied successfully to handle the demand- 
ing spatial resolution issue in nonlinear, time dependent DSA simulations. 
Berezhko and collaborators (e.g. ,[5]) developed a method that normalizes the 
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spatial variable by Xd{p) at each momentum value of interest during solu- 
tion of the CR kinetic equation. Their approach (which we here call 'nor- 
malized grid' or NG for convenience) allows efficient solution of the coupled 
system of time dependent gasdynamic equations and the CR transport equa- 
tion even when the diffusion coefficient has a strong momentum dependence 
(e.g.,K(p) oc p). It does depend on knowing the diffusion coefficient a priori 
of the simulations. The NG scheme was designed for simulations of supernova 
remnants, which have been represented by piston-driven spherical shocks in 
one-dimensional geometry. Assuming Bohm-like diffusion, Berezhko and El- 
lison [4,14] demonstrated qualitative agreement between nonlinear spherical 
DSA-modified shocks simulated with the above method and steady, plane 
Monte Carlo shocks that allow escape of CRs above a selected energy whose 
diffusion length corresponds to the curvature radius of the shock. On the other 
hand, since n{p) may evolve with the shock and may be nonuniform, one 
should have access to alternative methods that allow inclusion of non-steady 
diffusion properties and can still effectively handle the numerical challenges 
outlined above. 

As an alternative approach responding to these needs Kang et al. developed 
the CRASH (Cosmic-Ray Amr SHock) code in one dimensional (ID) plane- 
parallel geometry by combining Adaptive Mesh Refinement (AMR) techniques 
and a subgrid shock tracking technique [23] . The effectiveness of AMR in this 
situation comes from the fact that the highest resolutions are necessary only 
very close to the subshock, which can still be treated as a numerical discontinu- 
ity satisfying standard Rankine-Hugoniot relations, since the CR distribution 
is continuous across this feature. Consequently, individual refinement patches 
can be made quite small. By such efficient use of spatial gridding the blend 
of these computational strategies can greatly reduce the cost of time depen- 
dent DSA simulations, producing a practical tool for modeling CR modified 
shocks with arbitrary diffusive behavior. These benefits are especially valuable 
when combined with the so-called "Coarse-grained finite Momentum Volume" 
(CGMV) method for solving the diffusion-convection equation with a mini- 
mum of necessary information about the momentum dependence of the CR 
distribution function [18]. 

NG calculations of supernova remnants by Berezhko and collaborators (e.g. ,[5,6, 7,8]) 
have shown that, with Bohm-like diffusion, up to 50 % of explosion energy is 
converted to CRs when an assumed fraction 10~ 4 — 10~ 3 of incoming ther- 
mal particles are injected into the CR population at the subshock. Nonlin- 
ear modification to the flow structure can be substantial. Most obviously, as 
seen in many nonlinear CR shock simulations, the total density compression 
through the full shock structure is larger than the canonical gasdynamical 
value ("fg + l)/(7 fl — 1) = 4 for strong gasdynamical shocks. In SNR simu- 
lations by Berezhko and Volk [7,8], for example, 4 < P2/P0 < 15. Similar 
results obtain for nonlinear planar shocks. For example, we previously applied 
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our planar CRASH code to calculate the nonlinear evolution of CR modified 
shocks, also assuming Bohm-like CR diffusion [22,25]. Rather than assume a 
fixed injection fraction at the subshock, those simulations adopted a more so- 
phisticated "thermal leakage" injection model that filters subshock crossings by 
suprathermal particles using a momentum-dependent transparency function 
based on a nonlinear plasma model for postshock wave-particle interactions 
[27]. This injection model depends on a single, reasonably well-constrained 
parameter; namely, the ratio of the upstream magnetic field strength to down- 
stream turbulent magnetic field strength. Indeed in our simulations ~ 1CT 3 
of incoming thermal particles were injected into the CR population at strong 
quasi-parallel CR modified shocks. The ratio of CR energy to inflowing ki- 
netic energy increased with the shock Mach number, but approached ~ 0.5 
for large shock Mach numbers, M s > 30, and it was relatively independent of 
other upstream properties or variation in the injection parameter. The pres- 
ence of a preexisting, upstream CR population is equivalent to having slightly 
more efficient thermal leakage injection for such strong shocks, while it can 
substantially increase the overall CR energy in moderate strength shocks with 
M s < 3. 

In the current work we discuss an extension of the CRASH code to ID spher- 
ical geometry intended to model CR acceleration in SNR and spherical wind 
shocks. In order to implement shock tracking and AMR techniques effectively 
in this situation, we solve the fluid and diffusion-convection equations in a 
frame comoving with the outer spherical shock. Thus, the shock and refined 
region around it stay at the same grid location in the comoving frame, cre- 
ating a relatively straightforward scheme with all the advantages we found 
for plane geometry. The basic equations and details of the numerical method 
are described in §2. We will present simulation results for a representative 
Sedov- Taylor blast wave in §3, followed by a summary and comparison to 
previous simulation results in §4. Future papers will apply this code to the 
study of SNRs and their emissions as well as to quasi-spherical supersonic 
astrophysical winds. 



2 Numerical Method 

2.1 Basic Equations 

The evolution of CR modified shocks depends on a coupling between the 
gasdynamics and the CRs. That coupling takes place by way of resonant MHD 
waves, although it is common practice to express the pondermotive wave force 
and dissipation in the plasma using the associated CR pressure distribution 
properties along with a characteristic wave propagation speed (usually the 
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Alfven speed). Consequently, we solve the standard gasdynamic equations 
with CR pressure terms added in the conservative, Eulerian formulation for 
one dimensional spherically symmetric geometry. For strongly shocked flows 
numerical errors in computing the gas pressure from the total energy can lead 
to spurious entropy generation with standard methods, especially in the shock 
precursor. To avoid this problem we replace the usual energy conservation 
relation with the evolution of a modified entropy, S = Pg/p 13 ^ 1 , everywhere 
except across the subshock [24]. Energy conservation is applied across the 
subshock. The resulting dynamical equations are: 

dp 9 , , 2 . 

m + Tr ipu} = ~ pu ' (1) 

*gL + °^ + P, + P^-ltf, (2) 

+ ^ pegU + PgU) = ~ u l^~ ~ l( p6gU + PgU ^ ^ 

f + l(Su) = - 2 -Su + ^V(r,t) - L(r,t)}, (4) 



where P g and P c are the gas and the CR pressure, respectively, e g = -P g /[p(7 g — 1)] + 
u 2 /2 is the total energy of the gas per unit mass. The remaining variables, ex- 
cept for L and W have standard meanings. The injection energy loss term, 
L(r,t), accounts for the energy carried by the suprathermal particles injected 
into the CR component at the subshock and is subtracted from the post- 
shock gas immediately behind the subshock. Gas heating due to Alfv'en wave 
dissipation in the upstream region is represented by the term 

BP 

W(r,t) = -v A ^, (5) 

where va = B/y/Anp is the Alfv'en speed. This term derives from a simple 
model in which Alfv'en waves are amplified by streaming CRs and dissipated 
locally as heat in the precursor region (e.g. ; [31,17]). 

The CR population is evolved by solving the diffusion-convection equation , 
dg ( ,dg 1 d r 2/ , \y d 9 , n , 1 9 r 2 . dg 



where g = p 4 f, with f(p,r,t) the pitch angle averaged CR distribution, and 
where y = ln(p), while n(r,p) is the diffusion coefficient [34]. For simplicity 
we always express the particle momentum, p in units m p c and consider only 
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the proton CR component. The wave speed is set to be u w = va in the 
upstream region, while we use u w = in the downstream region. This term 
reflects the fact that the scattering by Alfv'en waves tends to isotropize the 
CR distribution in the wave frame rather than the gas frame [34]. Upstream, 
the waves are expected to be dominated by the streaming instability, so face 
upwind. Behind the shock various processes, including wave reflection, are 
expected to lead to a more nearly isotropic wave field (e.g., [2]). 

2.2 The CRASH code on a Comoving Spherical Grid 

Before introducing the spherical version of the CRASH code, it is helpful to 
outline the basic strategy of the analogous planar code. The full CR shock tran- 
sition involves a very wide range of length scales associated with the particle 
diffusion lengths that must be resolved to follow the shock evolution properly. 
In particular, structure must be resolved down close to the subshock physical 
thickness, since the diffusion scales for low energy, freshly injected particles are 
expected to be comparable. The outer scale of any DSA calculation is set by 
the greater of the full width of the shock transition or the size of the physical 
system being studied. Although it is necessary to include very fine scales to 
resolve interactions of the lowest energy CRs, their diffusion and acceleration 
are important over distances only within a few times Xd{pi n j) of the subshock, 
where Pi n j represents a characteristic injection momentum. Grid refinement 
in the CRASH code, therefore, includes a region spanning the subshock only 
large enough to include comfortably the diffusion scales of dynamically impor- 
tant CRs and enough levels to follow freshly injected CRs adequately [23]. In 
CRASH each level of refinement has twice the spatial resolution of the level 
below, but the same number of grid zones. Typically, we have used about 200 
zones in each refinement grid. To accomplish grid refinement effectively it is 
necessary to locate the subshock position exactly. Thus, we track the subshock 
as a moving, discontinuous jump inside the initial, uniform and fixed grid [23]. 
The subshock jump conditions and motion of the subshock are established by 
solution of the nonlinear gasdynamic Riemann problem. An additional mea- 
sure is used in the planar CRASH code to ensure that the shock remains near 
the middle of the refined region at all grid levels during each base-grid time 
step. In particular the fluid velocity in the refinement patch is transformed to 
the subshock rest frame calculated at the start of each base-grid time step. 
All the refined grids are then accordingly redefined at each time step. Thus, 
in the planar CRASH code the refinement patch moves smoothly with the 
shock, while the base grid is fixed in space. This greatly simplifies the refine- 
ment strategy, especially since it eliminates the need for construction of new 
refinement patches after the simulation is started. 

A somewhat different grid strategy is needed in spherical geometry, since the 



6 



analogous Galilean transformation cannot be applied. Drury and Mendonca 
[13] pointed out that a spherical shock can be made to be stationary by adopt- 
ing comoving variables which factor out a uniform expansion or contraction. 
Here we define a comoving frame that expands with the instantaneous shock 
speed. This strategy for simulations of SNR flows has, of course, been applied 
in various forms before (e.g., [10, 19]) and is the standard approach in cos- 
mological simulations (e.g., [33,30]). Following the conventional cosmological 
formalism, we adopt the comoving radial coordinate, x = r/a, where a is the 
expansion factor and a = 1 at the start of simulations. The expansion rate, 
a = {u s — v s )/x s , is found from the condition that the shock speed is zero 
at the comoving frame. Here u s and v s are the shock radial velocities in the 
Eulerian frame and in the comoving frame, respectively. Then the comoving 
density and pressures are defined as p = pa 3 , P g = P g a 3 , and P c = P c a 3 . The 
gasdynamic equation with CR pressure terms in the spherical comoving frame 
can be written as follows: 

dp Id 2 _ 

m + ad^ ipv) = -^ (7) 

^ + ~^-(pv 2 + P 9 + P e ) = ~-pv 2 - h -pv - dxp, (8) 
at a ox ax a 

d(pe Q ) 1 <9 . ~ . v dP c 2 ~ a _ . . 

— h - — {pe g v + P g v) = (pe g v + P g v) - 2-pe g - axpv, (9) 

dt adx a dx ax a 

I + \^ - ~> ~ # + i2 [W{x - t] ~ Ux - t)] - (10) 



As mentioned before, the entropy Eq. (10) is solved everywhere outside the 
subshock, while the energy conservation Eq. (9) is applied across the subshock. 
The deceleration rate is calculated numerically at a give time step by a = 
(a n — d n_1 )/At n , where d n_1 is the expansion rate at the previous time step. 
The diffusion-convection equation for the function g = p A f, is given by 



dt 



+ 



(v + u w ) 



1 d 

3ax 2 dx 



dx 



[x 2 (v + u w )} + -}( 



F -4 9 +3-5 + — -xki,?/-, 11 
ay a a 2 x 2 ox ox 



where f(p,r,t) is the comoving CR distribution function. This equation au- 
tomatically includes adiabatic losses for the CRs resulting from the spherical 
expansion of the flow. The equations 7-11 are solved using finite difference 
methods analogous to those outlined for the planar CRASH code and discussed 
more extensively in [23]. 
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2.3 The Thermal Leakage Injection Model 

In the "thermal leakage" model for CR injection at shocks, most of the down- 
stream thermal protons are locally confined by nonlinear MHD waves even 
when their speed exceeds the speed difference between the shock and the 
bulk downstream plasma. Only particles well into the tail of the postshock 
Maxwellian distribution can leak upstream across the subshock [27,28]. In 
particular, leaking particles not only must have velocities that would allow 
them to reach the receding subshock, they must also avoid being scattered by 
the MHD waves that mediate the plasma subshock. This filtering process is 
implemented numerically in the CRASH code by adopting a "transparency 
function", T esc (eB,v), that expresses the probability of supra-thermal particles 
at a given velocity, v , leaking upstream through the postshock MHD waves 
[16,24]. One free parameter controls this function; namely, e# = B Q /B^, which 
is the inverse ratio of the amplitude of the postshock MHD wave turbulence B± 
to the general magnetic field aligned with the shock normal, Bq [27]. Plasma 
hybrid simulations and theory both suggest that 0.25 < e# < 0.35 [27], so 
that the model is well constrained. We have also found for strong shocks that 
the time-asymptotic behaviors are only weakly dependent on the choice of 
this parameter. In this model the subshock is completely "opaque" to down- 
stream particles with momenta less than p±, (i.e., r esc = for p < p±, where 
Pi = (v>2/c)(1 + €b)/€b)- However, the subshock becomes virtually transparent 
to particles with p > (2 — 3)p 1 (i.e., r esc — > 1). We have ported this injection 
scheme to the spherical CRASH code, and will apply it to the simulations 
presented below. 



2.4 A Bohm-like Diffusion Model 

DSA, along with the time and length scales for CR acceleration, and subse- 
quent shock modification, depend on the spatial diffusion coefficient for the 
CRs. The diffusion coefficient can be expressed in terms of a mean scattering 
length, A, as n(x,p) = |At>. The Bohm diffusion model is commonly used to 
represent a saturated wave spectrum, providing the minimum expected diffu- 
sion coefficient in a quasi-parallel shock. We adopt a variation on that model 
in the simulations presented here. Bohm diffusion assumes that CRs scatter 
within one gyration radius (i.e., A ~ r g ). This gives k b oc p 2 /(p 2 + I) 1 / 2 with 
Kb oc p 2 in the limit of p « 1 and k b oc p in the limit of p » 1. Because the 
nonrelativistic momentum dependence is very steep, the diffusion coefficient 
and, thus, the diffusion length of freshly injected CRs, is very much smaller 
than for relativistic CRs if the shock speed is a small fraction of the speed of 
light. Since it is necessary, at least near the subshock, to resolve all relevant 
diffusion lengths, use of a Bohm diffusion model requires extremely fine grid 
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resolution where freshly injected CRs are concentrated. On the other hand, 
previous calculations have shown that, although details in diffusion models 
at the lowest momenta influence the thermal leakage injection rate and early 
shock evolution, they do not play significant roles once at least mildly relativis- 
ts particles dominate the CR pressure (e.g., [23]). Consequently, we will apply 
in the simulations presented here the following Bohm-like diffusion coefficient 
that includes a weaker nonrelativistic momentum dependence; namely, 



The coefficient n n = mc 3 /{3eB) = (3.13 x lO^cmV 1 )^- 1 , where B„ is the 
magnetic field strength in units of microgauss. The assumed density depen- 
dence for n accounts for compression of the perpendicular component of the 
wave magnetic field and also inhibits the acoustic instability in the precur- 
sor of highly modified CR shocks [21]. We note, also, that hereafter we use 
the subscripts '0', '1', and '2' to denote conditions far upstream of the shock, 
immediately upstream of the gas subshock and immediately downstream of 
the subshock, respectively. Thus, po represents the far-upstream gas density, 
which will be constant in the simulations presented here. 



3 Simulations of Sedov- Taylor Blast Waves 

For a supernova remnant (SNR) propagating into a uniform ISM (interstellar 
medium) the CR acceleration takes place mostly during free expansion and 
Sedov- Taylor (ST hereafter) stages, since the shock slows down significantly 
afterward. In fact the highest momentum, p max ~ 10 5 , which corresponds to 
the proton energy of 10 14 eV for a typical SNR, is roughly achieved by the 
end of the free expansion stage, while the transfer of explosion energy (E ) 
to the CR component occurs mostly during the ST stage. This is because the 
total CR energy gain is proportional to the kinetic energy passed through 
the shock, E sw = 2n f pou^r^dt, and E sw becomes comparable to E Q by the 
end of ST stage. In this paper, we calculate the CR acceleration during this 
stage, since the current version of CRASH code treats only the forward shock. 
Application of our new AMR algorithm for multiple spherical shocks is not 
simple, since it requires multiple, comoving grids. Although CR acceleration 
at the reverse shock can be dynamically important during the free expansion 
stage [15], CR energy generation occurs mostly at the forward shock during 
ST stage. Adiabatic losses by CRs accelerated early on in the interior and 
then advected outward through the ST phase would generally be very large. 




(12) 
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3. 1 Example SNR Model Parameters 



For most of the simulations we consider a supernova explosion with E Q = 
10 51 ergs and M sn = 10M Su n in a warm ISM with a uniform density, nn = 
0.3cm -3 and T = 10 4 K. The physical quantities are normalized, both in 
the numerical code and in the plots below, by the following constants: p Q = 
7.0 x 10~ 25 g cm' 3 , t = 1.3 x 10 3 yr, r Q = 6.14 pc, u Q = 4.6 x 10 3 km s~\ 
and P = 1.5 x 10~ 7 erg cm" 3 . We assume a Bohm-like diffusion coefficient, 
k = (6.26 x 10 21 cm 2 s" 1 )p(p /p), where P M = 5 is adopted for the mean 
ISM magnetic field. The Alfv'en speed is given by va = va,o/ \pI Po) 1 ^ 2 where 
va,o — 17 kms" 1 is for the background ISM. The pressure of the background 
gas is set to be P g o = 10~ 12 erg cm -3 so the gas sonic Mach number (M s = 
u s/ \JlgPgfll Po) of the initial shock is M s = 130. In order to explore effects of 
pre-existing CRs, we include an ambient (upstream) CR population, f{p) oc 
p" 4 - 5 , and set its pressure as either P C)0 = 0.005P 9>0 or P Cj0 = 0.5P Si0 - 111 
fact, the injection is so efficient for the strong shocks considered here that 
the presence of these pre-existing CRs is effectively equivalent to having an 
only slightly higher injection rate. It does not qualitatively affect the energy 
extracted by CRs from the shock or the degree of shock modification. 

Although the theoretically preferred values of e# lie between 0.25 and 0.35 
for strong shocks [26], values as large a these lead to very efficient initial 
thermal leakage injection for shock Mach numbers greater than 100 and a 
rate of shock evolution at simulation start-up that is very difficult to follow 
numerically. To avoid those start-up issues, we have to adopt a smaller value for 
the current simulations. Also, with the assumed k, thermal leakage injection of 
suprathermal particles is more efficient, compared to the true Bohm diffusion 
case. So we adopted e# ~ 0.16 in order to obtain a typical injection fraction of 
~ 10" 3 . Previously we showed that the CR acceleration depend only weakly 
on the values of e B [23], so this choice will not influence our conclusions. 

For the simulations presented here we use 230 uniformly spaced logarithmic 
momentum zones in the interval y = Inp = [lnpi,lnp2] when solving Eq. 
(11), where pi is determined instantaneously by the downstream flow speed, 
u 2 and lnj9 2 = 15.4. As the postshock density compression increases, u 2 and p\ 
decrease in time, so the particle distribution function g(y) from the previous 
time step solution is mapped onto the new momentum zones. 

The simulations are initialized at i = (t/t ) = 1 by the ST similarity solution, 
which is characterized by the shock position and speed expressed in code units, 
f s = (r s /r a ) = e s t 2/5 and u s = (u s /u ) = (2/5)^ 3/5 with & = 1.15167. The 
spatial resolution in the base grid is Afo = 6.0 x 10~ 4 , expressed in code units. 
In the 8th refined grid, which is the finest we have used in the simulations 
presented here, the resolution is Ar 8 = 2.3 x 10" 6 . Remarkably, we found 
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that although Afg is much larger than the diffusion length for p± 10 -2 , 
Idiff — 2.5 x 10 -9 , the spherical, comoving CRASH simulations show good 
numerical convergence. In particular, we found that an increase from 8 to 
10 levels of refined grids, corresponding to a factor 4 improvement in spatial 
resolution, leads to an increase of P c by less than 0.5 % (Convergence in P c 
is generally from below with our scheme.) Other variables respond similarly. 
This behavior is contrary to what we found in comparison simulations using a 
fixed, Eulerian grid. There, in agreement with our previous experience using 
Cartesian grids, spatial resolution several times finer than the diffusion length 
associated with p 1 were necessary for adequate convergence [23] . This will be 
discussed further below. 



3.2 Basic Code Tests 



As an initial test of the spherical CRASH code, we calculated identical ST 
blast waves on both a fixed Eulerian grid and a comoving grid. For the fixed 
Eulerian grid simulation our AMR scheme cannot be applied simply, because 
it requires that the refined patch move with the shock. Consequently, only 
the base grid (l g = 0) with Af = 6.0 x 10~ 4 is used in that simulation. 
Without AMR it is not practical to carry out that simulation with the high 
resolutions applied in subsequent AMR tests. We used instead the same base 
grid applied for the other simulations, but adopted a large diffusion coefficient, 
k = 0.1 p(po/p), to maintain adequate numerical convergence. In the comoving 
grid version of the test, two grid structures are considered; namely, the base 
grid alone, without AMR (i.e., the same as the fixed grid simulation except 
for the expansion), and the base grid with a single (l g = 1) refinement level 
using Afi = 3.0 x 10~ 4 . The large diffusion coefficient in this test model leads 
to very slow CR acceleration and associated shock modification. In order to 
introduce at least modest CR feedback effects the initial shock was much 
weaker than we used in follow-up tests, and the ambient CRs dominated the 
upstream pressure; namely, P Cj o = 4P Sj0 ? with f(p) oc p -4 ' 5 . The thermal 
leakage injection and Alfven wave terms were turned off for this test. We 
assumed a hot ISM in which n H = 3 x 10~ 3 cm -3 and T = 10 6 K with the 
initial shock Mach number M s = 13. The normalization constants are: p Q = 
7.0 x 10~ 27 g cur 3 , t = 6.1 x 10 3 yr, r Q = 28.5 pc, u Q = 4.6 x 10 3 km s -1 , and 
P Q = 1.5 x 10 _9 erg cm -3 . Fig. 1 shows that the three test simulations with 
different grid choices are in good agreement, demonstrating the validity of the 
comoving grid method. 

Our second test explored the convergence behavior of the spherical CRASH 
code. The physical problem was the same as for the first test, except that 
we used a much smaller and physically more interesting diffusion coefficient, 
k — 1.5 x 10~ 7 p(po/p) and a smaller ambient CR pressure, P Cj0 = 0.5P gfi . As 
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a benchmark, one simulation was carried out on a uniform, fixed grid iden- 
tical to that of the previous test. The CRASH, comoving simulations were 
done on a uniform grid both without AMR, l g = 0, and with AMR using 
l g = 1, 3, 5, 8, and 10. Results are shown in Fig. 2. The AMR simulation 
results converge well as the number of refinement levels increases. CR ac- 
celeration efficiency increases with finer resolution, consistent with previous 
results (e.g., [23]). As expected, the coarse fixed grid simulation and the co- 
moving simulations with coarser resolution are not very accurate as measured 
by the high resolution solutions. It is interesting to note, however, that the 
CR acceleration is more accurate in the CRASH simulation with the comov- 
ing grid (with l g = 0) than the fixed Eulerian grid simulation of comparable 
resolution. Thus, the comoving grid provides a more economical solution to 
the problem than the fixed grid. This can be understood as follows. In the 
comoving grid simulation, the peak of the spatial distribution of f(p) and the 
subshock location remain in the same grid zone throughout the calculation. 
Consequently, the subshock compression (V • u), which determines the rate at 
which low energy CRs accelerate (see Eq. [11]), is applied to the CR distribu- 
tion consistently. In the fixed Eulerian grid, on the other hand, the shock drifts 
through the grid as the shock evolves. As this happens the numerical estimate 
of subshock compression varies and, on average is reduced. This reduces the 
effective adiabatic compression applied to f(p) at the subshock. This behavior 
also leads to the the unsteady behavior of the CR pressure at the subshock, 
P Cj 2, for the fixed grid simulation, shown in the lower right panel of Fig. 2. 



3.3 Results 



Figs. 3 and 4 show the early evolution of the ST blast wave outlined in §3.1 
(initial Mach number 130) over the interval 1.0 < t < 1.5. The simulation 
shown in Fig. 3 includes the Alfven wave terms, while the simulation shown 
in Fig. 4 turns them off. Otherwise they are identical. The black solid lines 
show the initial similarity solution with P c0 = 0.5P 9i o at f = 1. The CR 
pressure reaches the highest value during the initial 130 years (At = 0.1), 
but it decreases as the shock slows down over time . As previously shown 
in SNR simulations by Berezhko et al. [7], and anticipated from earlier work 
(e.g., [31,1,29,17]), the simulation with Alfv'en wave drift and heating terms 
(Fig. 3) shows less efficient CR acceleration than the simulation without those 
terms (Fig. 4). After 650 years (At = 0.5), p max ~ 10 4 , as indicated by 
the exponential cut-off in the particle distribution function g(p) at the shock 
shown in the lower left panel of Figs. 3 and 4. In both simulations, nonlinear 
modifications due to CRs are significant, represented clearly by growth of a 
significant precursor and concave curvatures in the CR particle distribution. 

The density in the precursor just in front of the subshock, p 1 , and in the post- 
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shock region, p2, immediately after the subshock, respectively, are shown in 
the top panels of Fig. 5. In the model with the Alfven wave terms, the pre- 
cursor compression stays around 2, while the total compression stays around 
7-8. Without the wave terms, on the other hand, the CR acceleration is very 
efficient and the total density compression can go up to a factor of 30. This last 
result is consistent with our simulations of plane-parallel shocks omitting wave 
terms, i.e., P2/P0 ~ 1.5M s °q where M Si o is the Mach number of a unmodified 
shock [24] as well as previous SNR simulations without these terms (e.g., [6]). 
The second from the top panels show the Mach number of the subshock, which 
decreases quickly during the initial phase and stays around M s pa 5 for the 
model with the wave terms. The third from the top panels show the the CR 
pressure and gas pressure normalized to the ram pressure of the unmodified 
ST similarity solution, i.e., PqU 2 st oc f~ 6 / 5 . After the initial quick changes, 
these ratios remain constant at P Ci2 /(po^|t) ~ 0-4 an d -P s ,2/(Po m |t) ~ 0-2 
in the model with the wave terms. Consequently, the conversion factor from 
the shock kinetic energy to the CR energy approaches an asymptotic value, 
although the pressure itself decreases as the remnant expands. The bottom 
panels show the fraction of particles that have passed through the shock and 
injected into the CR component. With the adopted value of e# = 0.16, this 
fraction is £ ~ 2 x 10~ 4 — 2 x 10" 3 in the model with the wave terms, while 
it can be quite large, £ ~ 0.01, during the initial phase in the model without 
the wave terms. 

Fig. 6 shows the total CR numbers integrated over the simulated volume, 
N P = J r 2 drf(p)p 2 , and G p = N p p 2 = J r 2 drg(p) in code units. The slope of 
the integrated spectrum, q = —d(\nN p )/d\np, is shown in the bottom panels. 
We note that it is much easier to recognize the nonlinear concave curvature in 
the plot of G p , while N p resembles a power law below p m ax- In the model with 
wave terms, the slope for nonrelativistic particles is q ~ 2.2, consistent with 
the compression ratio across the subshock with M s pa 4.6. The slope flattens 
to q ~ 1.6 near the highest momentum, above which the spectrum decreases 
exponentially. This agrees well with Berezhko and Volk [7] where the power 
law index of N p is about 1.7 for p > 10 3 , and steeper at lower energy for 
the model with their injection rate, 77 = 10~ 3 . The steep momentum slope 
at low energies results in nonlinear CR shocks from the fact that these CRs 
are accelerated primarily through localized crossings of the subshock. The flat 
high energy slope, on the other hand, reflects the fact that these CRs, with 
their large mean free scattering paths, are accelerated across the full shock 
precursor, so the relevant compression is greater. 

Finally, we show in Fig. 7 the total kinetic, thermal and CR energies integrated 
over the simulation volume, E^ n , E t h, Ec, respectively. This shows that the 
fraction of total kinetic energy remains the same as that of the initial similarity 
solution. But, by t = 10 the gas thermal energy fraction decreases as the CR 
energy fraction increase to ~ 60 %. 
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The SNR models for our calculations are very similar to those presented by 
Berezhko and Volk [7] except that they started their calculations from the 
free expansion stage, and their supernova ejecta dynamics was represented by 
a piston using the so-called thin-shell approximation. Once the ST stage is 
underway, these different initial conditions are not very important to the CR 
evolution. But their numerical technique is quite different from ours, as men- 
tioned in §1. The results of their model with the fixed injection rate rj = 10~ 3 
and Alfven wave damping is similar to those of our model with Alfven wave 
terms included. The overall results, that is, density compression both across 
the subshock and across the entire shock, the integrated particle spectrum, 
and the total CR energy fraction are all in good agreement. As expected, the 
early evolution of our model is somewhat different, because we started at the 
beginning of the ST stage (t = 1), while they began with the earlier free 
expansion stage (t — 0). 



4 Summary 

We have developed a new cosmic ray shock code in one dimensional spherical 
geometry. The diffusion convection equation for the CR distribution function is 
solved along with the gas dynamics equations modified to include the diffusive 
CR pressure. Simple models for thermal leakage injection and Alfv'en wave 
propagation and dissipation are also included. In order to implement the shock 
tracking and AMR techniques, we adopt a comoving spherical grid which 
expands with the instantaneous shock speed. In the comoving grid, the shock 
remains at the same location, so the compression rate is applied consistently 
to the CR distribution at the subshock. This results in much more accurate 
and efficient low energy CR acceleration and faster numerical convergence on 
coarser grid spacings, compared to the simulations in a fixed, Eulerian grid. 

We have calculated the CR acceleration during the Sedov- Taylor stage for a 
typical supernova remnant in a warm interstellar medium. The Mach number 
of the initial shock is M s = 130. In the simulations without the Alfven wave 
drift and dissipation terms, the results are consistent with what we found 
for shocks with similar Mach numbers in plane-parallel geometry reported in 
previous papers [24,25]. Since the diffusion length of the highest momentum 
(Pmax/ m c ~ 10 6 ), Idiff / r o ~ 0.25, is still smaller than the shock radius, r s /r Q ~ 
2.9, at t/t = 10, the pressure dilution in the spherically expanding volume 
should have only minor effects. Although the shock expands and slows down, 
the CR acceleration efficiency expressed in terms of P c ,2/[poUst\ * s similar to 
that for plane-parallel shocks with similar initial shock parameters. 

When the drift and dissipation of Alfv'en waves in the upstream region are 
included, the CR acceleration is less efficient, as expected from previous studies 
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[31,29,17,7]. We summarize here the main properties of the simulations with 
such wave effects included: 

1. The postshock pressures, P c and P g , normalized by the unmodified ST shock 
ram pressure, PoU 2 ST oc £~ 6 / 5 , approach time- asymptotic values quickly. The 
postshock CR pressure is about 40 % of the ST shock ram pressure, while the 
gas pressure drops from its gasdynamic value of 75 % to about 20 % of the 
shock ram pressure. 

2. A significant shock precursor develops in response to nonlinear feedback 
from the CR pressure, so the subshock weakens to M s 5. The density 
compression factor through the precursor is pi/po ~ 2, while the compression 
over the total transition is P2/P0 ~ 8- However, as seen in previous spherical 
CR shock studies (e.g., [12,7]), the subshock does not completely disappear, as 
in the case of plane-parallel shock simulations with the self-adjusting thermal 
leakage injection model [24]. 

3. Both the CR momentum distribution at the shock, g(r s ,p) = f(r s ,p)p i , and 
the integrated distribution, G(p) = J r 2 drg(p), clearly exhibit characteristic 
concave curvature, reflecting the nonlinear velocity structure in the precursor. 
On the other hand, the integrated number distribution, N p (p) = J r 2 drf(p)p 2 
does not reveal the nonlinear features so obviously. At nonrelativistic momenta 
N p oc p~ 2 ' 2 , while it flattens to N p oc p~ ls at the highest momenta. 

As a basic test of our new scheme for studying spherical, CR modified shock 
evolution, we simulated SNR blast waves that were similar to SNR models 
published by Berezhko & Volk [7] based on very different numerical methods. 
The results are in substantial agreement. Our new CRASH code in a comov- 
ing spherical code made this comparison possible for the first time, because 
the calculation of DSA in a comoving spherical grid achieves numerical con- 
vergence at a grid resolution much coarser than that required in an Eulerian 
grid. 
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Fig. 1. Numerical convergence test for a Sedov- Taylor blast wave in a hot ISM 
with uh = 3 x 10~ 3 cm~ 3 and T = 10 6 K. The initial Mach number of the 
shock is M s = 13. A preexisting CR population, f\p) oc p~ 4 ' 5 , corresponding to 
P C: Q = 4P Si o> is assumed. The diffusion model is k{p) = 0.1p(po/p) in code units. 
The black dashed lines show the initial conditions. The CR spectrum at the shock , 
g(r s ,p) = f(r s ,p)p 4: , is shown in the lower left panel, while the postshock CR pres- 
sure normalized by the fiducial pressure is shown as a function of time in the lower 
right panel. Black lines are for the calculation in the fixed Eulerian grid without 
AMR, blue lines for the calculation in the comoving grid without AMR, and red 
lines for the calculation in the comoving grid with one level of AMR. 
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Fig. 2. Numerical convergence test for a Sedov- Taylor blast wave in a hot ISM with 
fin = 3 x 10~ 3 cm~ 3 and T = 10 6 K. A preexisting CR population, f(p) oc p -4,5 , 
corresponding to P Cj o = 0.5-P 9 .o, is assumed (yellow dotted line in the low right 
panel). The diffusion model is R{p) = 1.5 x W~ 7 p(po/p). The CR spectrum at the 
shock, g(r s ,p) = f(r s ,p)p 4 , is shown in the lower left panel, while the postshock 
CR pressure is shown as a function of time in the lower right panel. Black dotted 
lines are used for the calculation with the fixed Eulerian grid without AMR. For 
the calculations with the comoving grid, different colored lines are used for different 
levels of AMR: cyan, magenta, green, blue, red and black for l g = 0, 1, 3, 5, 8, and 
10, respectively. 
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Fig. 3. Evolution a SNR expanding into the warm uniform ISM with nn = 0.3cm -3 
and T = 10 4 K. The model parameters are E Q = 10 51 ergs, M sn = 10Ms U n, and 
= 5. The model assumes a preexisting CR population of f(p) oc p~ 4 5 , with 
PcO = 0.5P 9i o- The injection parameter for thermal leakage injection, €b = 0.16. 
The lower left panel shows the CR spectrum at the shock, g(r s ,p) = f(r s ,p)p 4 . The 
initial condition at t = 1.0 (solid line) is set by the Sedov- Taylor similarity solution. 
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